Method and System for Modeling Fractures in Ductile Rock

ABSTRACT

Method and systems for modeling fractures in quasi-brittle materials are provided. An exemplary method included generating a model that incorporates a unified creep-plasticity (UCP) representation into a constitutive model for a ductile rock. The model may be used in a finite element analysis to model hydraulic fractures in the ductile rock.

This application claims the benefit of U.S. Provisional PatentApplication 61/359,157 filed Jun. 28, 2010 entitled METHOD AND SYSTEMFOR FRACTURES IN DUCTILE ROCK, the entirety of which is incorporated byreference herein.

FIELD OF THE INVENTION

Exemplary embodiments of the present techniques relate to a method andsystem for the analysis and modeling of fractures in quasi-brittlematerials. Specifically, an exemplary embodiment provides for modelinghydraulic fractures in ductile rock.

BACKGROUND

This section is intended to introduce various aspects of the art, whichmay be associated with exemplary embodiments of the present techniques.This discussion is believed to assist in providing a framework tofacilitate a better understanding of particular aspects of the presenttechniques. Accordingly, it should be understood that this sectionshould be read in this light, and not necessarily as admissions of priorart.

As hydrocarbon reservoirs that are easily harvested, such as oilreservoirs on land or reservoirs located in shallow ocean water, areused up, other hydrocarbon sources must be used to keep up with energydemands. Such reservoirs may include any number of unconventionalhydrocarbon sources, such as biomass, deep-water oil reservoirs, andnatural gas from other sources.

One such unconventional hydrocarbon source is natural gas produced fromshale, termed “shale gas.” Because shale may have insufficientpermeability to allow significant fluid flow to a well bore, many shalesare currently not considered as commercial sources of natural gas.However, shale gas has been produced for years from shales with naturalfractures. Recently, a significant increase in shale gas production hasresulted from hydraulic fracturing, which is used to create extensiveartificial fractures around well bores. When combined with horizontaldrilling, which is often used with shale gas wells, the hydraulicfracturing may allow formerly unpractical shale layers to becommercially viable.

Shales that host economic quantities of gas have a number of commonproperties. They are rich in organic material and are usually maturepetroleum source rocks in the thermogenic gas window. They are oftensufficiently brittle and rigid enough to maintain open fractures, over ashort period. Some of the gas produced is held in natural fractures,some in pore spaces, and some is adsorbed onto the organic material. Thegas in the fractures can be produced immediately, while the gas adsorbedonto organic materials may be released as the formation pressuredeclines.

The fracturing process is complicated and can be modeled to improve theefficiency of the fracturing. A more efficient fracturing process maylead to a more productive reservoir. In other words, a greater amount ofthe gas trapped in a shale layer may be harvested. As discussed below, anumber of researchers have proposed different techniques for modelingthe growth of fractures in various rock formations.

In Perkins, T. K., and Kern, L. R., “Widths of Hydraulic Fractures,”JPT, September, 222 Trans. AIME 937-949 (1961), a fluid mechanics studywas performed on the rupture of brittle materials and the theory ofelastic deformation of rock. The results indicated that for a givenformation, crack (or fracture) width is essentially controlled by fluidpressure drop. High pressure drops result in relatively wide cracks,while low pressure drops result in relatively narrow cracks.

In Rice, J. R. and Rosengren, G. F., “Plane strain deformation near acrack tip in a power-law hardening material,” 16 J. Mech. Phys. Solids1-12 (1968), crack-tip strain singularities were investigated with theaid of an energy line integral exhibiting path independence for allcontours surrounding a crack tip in a two-dimensional deformation field.Elastic and elastic/plastic materials were studied. The study determinedthat the product of stress and strain exhibits a singularity varyinginversely with distance from the tip in all materials.

In Nordgren, R. P., “A Propagation of a Vertical Hydraulic Fracture,”SPEJ, August, 306-314, Trans., AIME, 253 (1972), the propagation ofhydraulic fractures of limited vertical and elliptic cross-section werestudied with the inclusion of fluid loss. The study determined that thefracture length and width grow faster with time in the no-loss than inthe large loss case.

In Leung, C. K. Y., and Li, V. C., 24 Journal of Materials Science854-862 (1989), an experimental technique for indirectly determining atension-softening curve, was developed. The tension-softening curve issuggested to be a size-independent fracture ‘parameter’ forquasi-brittle materials. The technique was based on generating specimensthat had a large size and a four-point bending mode.

In McDowell D. L., Miller M. P., and Brooks D.C., “A UnifiedCreep-Plasticity Theory for Solder Alloys”, 42 Fatigue ElectronicMaterial AST STP 1153 (1994), the fracture behavior of solders wasinvestigated. The goals of the study were to develop a thermo-mechanicalunified creep-plasticity model, which could be used to characterize theresponse of solders, and to characterize a selected solder alloy throughmechanical testing and correlation with the model across a range ofstrain rates, temperatures, aging times, and loading-unloadingconditions. The model allowed the performance of solders under fatigueand creep-plasticity interactions due to thermo-mechanical cycling to bemodeled.

In Benzeggagh and Kenane, “Measurement of mixed-mode delaminationfracture toughness of unidirectional glass/epoxy composites withmixed-mode bending apparatus,” 56 Composites Sci. Tech. 439-449 (1996),the initiation of cracking and delamination growth in a unidirectionalglass/epoxy composite were evaluated under mode I, mode II, and mixedmode I+II static loading. In addition to the characterization ofdelamination, the apparatus allowed the plotting of an R curve, fromwhich a total fracture resistance, G_(TR), could be obtained.

In Martin, A. N., “Crack Tip Plasticity: A Different Approach toModelling Fracture Propagation in Soft Formations,” SPE 63171 (2000), acrack tip plasticity (CTP) method was developed, which assumed afracture tip of finite radius, with a zone of plastically deformedmaterial around it. The plastic zone forms when stresses on the rockincrease beyond the yield point of the rock, at which point the rockstarts to flow plastically. The plastic zone acts to absorb extra energyfrom the fracturing fluid, making it harder to propagate fracturesthrough formations with significant plastic deformation. This in turnmeans that, for a given ductile material, fractures will be smaller andless conductive than those predicted by LEFM.

In van Dam, D. B., “Impact of Rock Plasticity on Hydraulic FracturePropagation and Closure,” SPE 63172 (2000), scaled laboratoryexperiments of hydraulic fracture propagation and closure in softartificial rock and outcrop rock samples were performed. Numericalsimulations of the fracture behavior in plastic rocks were alsoperformed, using independently measured rock properties. The simulationsaided in interpreting the measurements and extrapolating the results tofield scale. Plasticity induces a larger width in a fracture for a givennet pressure, compared with elastic rock. However, the pressure topropagate fractures was only marginally increased and, in the case oflaboratory tests, was actually lower than expected from elasticbehavior. The most dramatic effect of plasticity is that closure is muchlower than the confining stress due to strong stress redistributionalong the fracture.

In Dean, R. H. and Schmidt, J. H., “Hydraulic Fracture Predictions witha Fully Coupled Geomechanical Reservoir Simulator,” SPE 116470 (2008), ageomechanical reservoir simulator was developed that combined hydraulicfracture growth, multiphase/multicomponent Darcy/nonDarcy porous flow,heat convection and conduction, solids deposition, andporoelastic/poroplastic deformation in a single application. The programcontained two separate criteria that could be used to model fracturepropagation: a critical fracture-opening criterion based on a stressintensity factor and a cohesive zone model that used quadrilateralcohesive elements in the fracture. The cohesive zone model includes acohesive strength and an energy release rate in the calculations at thetip of a propagating hydraulic fracture.

U.S. Patent Application Publication No. 2008/0091396 by Kennon, et al.,provides a method and system for modeling and predicting hydraulicfracture performance in hydrocarbon reservoirs. An unstructuredautomatic mesh is generated and computational algorithms are executedusing a finite element numerical approach. The method is to model ahydrocarbon reservoir, wells, and completions as a single system,accounting for static information and transient behavior of wells,hydraulic fractures, and reservoirs in a single model. Models used forthe fracture behavior are not disclosed.

Most of the current hydraulic fracture modeling in the oil industrycontinues to rely on empirical models, or on numerical methods based onlinear elastic fracture mechanics (LEFM). Generally, these methods arevalid for brittle materials and give reasonable predictions forhydraulic fractures in hard rocks. However, for ductile rock, such asshale, clay, weakly consolidated sandstones (low cohesion granularmaterials), or other relatively less brittle, softer, and/or moredeformable rocks, the LEFM-based methods typically give overlyconservative and often limited-use or unuseful predictions of fracturegeometry.

SUMMARY

For ductile rock (which are quasi-brittle in nature), cohesive fracturemechanics based method can be used to simulate the fracture tip effect.

An exemplary embodiment of the present techniques provides a method formodeling fractures. The method includes generating a model thatincorporates a unified creep-plasticity (UCP) representation into aconstitutive model for a ductile rock.

The ductile rock forms fractures that may have a fracture process zonedistributed over a finite size volume, which exhibits strain-softening.A ratio of a length of the fracture process zone to a length of anon-linear fracture zone located at the tip of the fracture may bebetween about 1:1 and 1:1.5.

The method may include predicting a J-integral value and stress fieldincorporated plasticity effect from finite element analysis using themodel. The method may further include determining the effective fractureproperties of the ductile rock. Further, the method may includedetermining the effective fracture properties of the ductile rock usinga numerical method. A fracture in the ductile rock may be predicted bythe method. A hydraulic fracture of the ductile rock may be predictedusing a pore pressure cohesive zone model.

A flow of fluid may be determined in a fractured ductile rock using themethod. A parameter in a fracturing process may be adjusted based, atleast in part, on the predicted fracturing of the ductile rock. Theparameter may be a mass flow rate, a pressure, or a combination thereof.

Another exemplary embodiment provides a method for modeling hydraulicfractures in ductile rocks. The method includes generating aconstitutive model for a quasi-brittle material that incorporatesunified creep-plasticity, performing a finite element analysis using theconstitutive model, determining an effective fracture property for thequasi-brittle material, and predicting hydraulic fractures in thequasi-brittle material.

The method may also include adjusting a hydraulic fracturing parameterbased, at least in part, on the fractures predicted in the ductile rock.A ratio between a length of a fracture process zone and a length of anon-linear fracture zone at the tip of a fracture may be between about1:1 and 1:1.5.

Another exemplary embodiment provides a system for modeling fractures ina reservoir. The system includes a processor and a storage mediumcomprising a data structure comprising a model that incorporates unifiedcreep-plasticity (UCP) into a constitutive model for a ductile rock. Thesystem also includes a computer readable medium that includes codeconfigured to direct the processor to perform a finite element analysisof the ductile rock using the model.

The processor may include a multi-processor cluster computer system. Thecomputer readable medium may include code configured to direct theprocessor to run a finite element analysis for predicting a J-integralvalue and a stress field incorporating plasticity effect using themodel. The computer readable medium may also include code configured todirect the processor to determine the fracture properties of the ductilerock. The computer readable medium may include code configured to directthe processor to predict a hydraulic fracture in the ductile rock.Further, the computer readable medium may include code configured todirect the processor to determine a flow of hydrocarbon in a fracturedductile rock.

Another exemplary embodiment of the present techniques provides anon-transitory, computer-readable medium that includes a data structurerepresenting a model that combines a constitutive model with a unifiedcreep-plasticity model. The non-transitory, computer-readable mayinclude code configured to direct a processor to perform a finiteelement analysis using the model, determine an effective fractureproperty for a ductile rock, and predict fractures in the ductile rock.

Another exemplary embodiment provides a method for harvestinghydrocarbons from a ductile rock formation. The method includesgenerating a model that incorporates unified creep-plasticity (UCP) intoa constitutive model for a ductile rock, performing a finite elementanalysis using the model, and predicting fractures in a ductile rockformation based, at least in part, upon a result from the finite elementanalysis. The method also includes fracturing the ductile rock formationusing parameters based, at least in part, on the prediction of thefractures. Further, the method includes harvesting hydrocarbons from thefractured rock.

BRIEF DESCRIPTION OF THE DRAWINGS

The advantages of the present techniques are better understood byreferring to the following detailed description and the attacheddrawings, in which:

FIG. 1 is a diagram of a hydrocarbon reservoir, in accordance with anexemplary embodiment of the present techniques;

FIG. 2 is a diagram that shows hydraulic fracturing through a singlewell extending from a surface facility through an overburden into areservoir, in accordance with an exemplary embodiment of the presenttechniques;

FIG. 3 is a close-up view of a fracture forming in the reservoir of FIG.2;

FIG. 4 shows a non-linear fracture zone ahead of a fracture in materialsranging from brittle to ductile;

FIG. 5 is a process flow chart of a method for modeling fractures inductile rock, in accordance with an exemplary embodiment of the presenttechniques;

FIG. 6 is a graph of a von Mises stress distribution ahead of a fracturetip, in accordance with an exemplary embodiment of the presenttechniques;

FIG. 7 is graph illustrating the estimation of the fracture toughnessfor a nonlinear softening material-ductile rock based on the size of thefracture process zone, in accordance with an exemplary embodiment of thepresent techniques;

FIG. 8 is a graph illustrating the cohesive element traction-separation(T-S) Law used in the CZM, in accordance with an exemplary embodiment ofthe present technique;

FIG. 9 is a drawing of a three-layer 3D pore pressure cohesive elementmodel that may be used to predict hydraulic fracture caused by fluidinjection using exemplary embodiments of the current techniques;

FIG. 10 is an illustration of a fracture geometry predicted for themodel of FIG. 9, in accordance with exemplary embodiments of the presenttechniques;

FIG. 11 is a series of graphs illustrating predicted fracturepropagation prediction of finite element analysis using a UCPconstitutive model, in accordance with exemplary embodiments of thepresent techniques;

FIG. 12 is a series of graphs from a finite element analysis thatillustrate changes in the fracture propagation predictions from a pseudo3D model, in accordance with exemplary embodiments of the presenttechniques;

FIG. 13 is a graph of effective fracture toughness at different leak-offcoefficients, in accordance with exemplary embodiments of the presenttechniques; and

FIG. 14 is a block diagram of an exemplary cluster computing system 1400that may be used to implement the present techniques, in accordance withexemplary embodiments of the present techniques.

DETAILED DESCRIPTION

In the following detailed description section, the specific embodimentsof the present techniques are described in connection with preferredembodiments. However, to the extent that the following description isspecific to a particular embodiment or a particular use of the presenttechniques, this is intended to be for exemplary purposes only andsimply provides a description of the exemplary embodiments. Accordingly,the present techniques are not limited to the specific embodimentsdescribed below, but rather, such techniques include all alternatives,modifications, and equivalents falling within the true spirit and scopeof the appended claims.

At the outset, and for ease of reference, certain terms used in thisapplication and their meanings as used in this context are set forth. Tothe extent a term used herein is not defined below, it should be giventhe broadest definition persons in the pertinent art have given thatterm as reflected in at least one printed publication or issued patent.Further, the present techniques are not limited by the usage of theterms shown below, as all equivalents, synonyms, new developments, andterms or techniques that serve the same or a similar purpose areconsidered to be within the scope of the present claims.

“Algorithm” carries its normal meaning and refers without limitation toany series of repeatable steps that result in a discrete value orvalues. For example, an algorithm may include any mathematical,statistical, positional, or relational calculation between any numbersof user-specified, preset, automatically-determined, or industry- orsystem-acceptable data elements. In several embodiments, variousalgorithms may be performed on subject data elements in relation to apreviously defined data evaluation sample in order to produce a single,meaningful data value.

“Computer-readable medium” or “Non-transitory, computer-readable medium”as used herein refers to any non-transitory storage medium thatparticipates in providing instructions to a processor for execution.Such a medium may include both non-volatile media and volatile media.Non-volatile media includes, for example, NVRAM, or magnetic or opticaldisks. Volatile media includes dynamic memory, such as main memory.Common forms of non-transitory, computer-readable media include, forexample, a floppy disk, a flexible disk, a hard disk, an array of harddisks, a magnetic tape, or any other magnetic medium, magneto-opticalmedium, a CD-ROM, a holographic medium, and any other optical mediumfrom which a computer can read data or instructions. Further, thenon-transitory, computer-readable media may also include a memorydevice, such as a RAM, a PROM, an EPROM, a FLASH-EPROM, a solid statemedium like a memory card, any other memory chip or cartridge, or anyother tangible medium from which a computer can read data orinstructions.

A “rock constitutive model” is a data construct that interrelates stressand strain properties. In physics, a constitutive equation is a relationbetween two physical quantities (often described by tensors) that isspecific to a material or substance, and approximates the response ofthat material to external forces. It is combined with other equationsgoverning physical laws to solve physical problems. For the analysisdescribed herein, constitutive relations connect applied stresses orforces to strains or deformations. The stress-strain constitutiverelation for linear materials commonly is known as Hooke's law. Inexemplary embodiments, the rock constitutive model described hereinincorporates nonlinear effects, strain rate, and temperaturedependencies, as well as plasticity and creep effects. The model may bestored in a non-transitory computer-readable medium for use by aprocessor during simulation and modeling calculations.

“Exemplary” is used exclusively herein to mean “serving as an example,instance, or illustration.” Any embodiment described herein as“exemplary” is not to be construed as preferred or advantageous overother embodiments.

“Formation” refers to a body of rock or other subsurface solids that issufficiently distinctive and continuous that it can be mapped, forexample, by seismic techniques. A formation can be a body of rock ofpredominantly one type or a combination of types. A formation cancontain one or more hydrocarbon-bearing zones. Note that the termsformation, reservoir, and interval may be used interchangeably, but willgenerally be used to denote progressively smaller subsurface regions,zones, or volumes. More specifically, a formation will generally be thelargest subsurface region, a reservoir will generally be a region withinthe formation and will generally be a hydrocarbon-bearing zone (aformation, reservoir, or interval having oil, gas, heavy oil, and anycombination thereof), and an interval will generally refer to asub-region or portion of a reservoir. A hydrocarbon-bearing zone can beseparated from other hydrocarbon-bearing zones by zones of lowerpermeability such as mudstones, shales, or shale-like (highly compacted)sands. In one or more embodiments, a hydrocarbon-bearing zone includesheavy oil in addition to sand, clay, or other porous solids.

A “fracture” is a crack or surface of breakage within rock not relatedto foliation or cleavage in metamorphic rock along which there has beenminimal movement. A fracture along which there has been lateraldisplacement may be termed a fault. When walls of a fracture have movedonly normal to each other, the fracture may be termed a joint. Fracturesmay enhance permeability of rocks greatly by connecting pores together,and for that reason, joints and faults may be induced mechanically insome reservoirs in order to increase fluid flow.

“Fracturing” refers to the structural degradation of a subsurface shaleformation from applied thermal or mechanical stress. Such structuraldegradation generally enhances the permeability of the shale to fluidsand increases the accessibility of the hydrocarbon component to suchfluids.

“Fracture gradient” refers to a fluid pressure sufficient to create orenhance one or more fractures in the subterranean formation. As usedherein, the “fracture gradient” of a layered formation also encompassesa parting fluid pressure sufficient to separate two or more adjacentbedding planes in a layered formation. It should be understood that aperson of ordinary skill in the art could perform a simple leak-off teston a core sample of a formation to determine the fracture gradient of aparticular formation.

“Hydrocarbons” are generally defined as molecules formed primarily ofcarbon and hydrogen atoms such as oil and natural gas. Hydrocarbons mayalso include other elements, such as, but not limited to, halogens,metallic elements, nitrogen, oxygen, and/or sulfur. Hydrocarbons may beproduced from hydrocarbon reservoirs through wells penetrating ahydrocarbon containing formation. Hydrocarbons derived from ahydrocarbon reservoir may include, but are not limited to, kerogen,bitumen, pyrobitumen, asphaltenes, oils, or combinations thereof.Hydrocarbons may be located within or adjacent to mineral matriceswithin the earth. Matrices may include, but are not limited to,sedimentary rock, sands, silicilytes, carbonates, diatomites, and otherporous media.

A “hydraulic fracture” is a fracture at least partially propagated intoa formation, wherein the fracture is created through injection ofpressurized fluids into the formation. While the term “hydraulicfracture” is used, the techniques described herein are not limited touse in hydraulic fractures. The techniques may be suitable for use inany fracture created in any manner considered suitable by one skilled inthe art. Hydraulic fractures may be substantially horizontal inorientation, substantially vertical in orientation, or oriented alongany other plane. Generally, the fractures tend to be vertical at greaterdepths, due to the increased pressure of the overburden.

“Hydraulic fracturing” is a process used to create fractures that extendfrom the well bore into formations to stimulate the potential forproduction. A fracturing fluid, typically viscous, is generally injectedinto the formation with sufficient pressure (for example, at a pressuregreater than the lithostatic pressure of the formation) to create andextend a fracture. A proppant may often be used to “prop” or hold openthe created fracture after the hydraulic pressure used to generate thefracture has been released. Parameters that may be useful forcontrolling the fracturing process include the pressure of the hydraulicfluid, the viscosity of the hydraulic fluid, the mass flow rate of thehydraulic fluid, the amount of proppant, and the like.

As used herein, “material properties” represents any number of physicalconstants that reflect the behavior of a ductile rock or quasi-brittlematerial. Such material properties may include, for example, Young'smodulus (E), Poisson's Ratio (Φ), tensile strength, compressionstrength, shear strength, creep behavior, and other properties. Thematerial properties may be measured by any combinations of tests,including, among others, a “Standard Test Method for UnconfinedCompressive Strength of Intact Rock Core Specimens,” ASTM D 2938-95; a“Standard Test Method for Splitting Tensile Strength of Intact Rock CoreSpecimens [Brazilian Method],” ASTM D 3967-95a Reapproved 1992; a“Standard Test Method for Determination of the Point Load Strength Indexof Rock,” ASTM D 5731-95; “Standard Practices for Preparing Rock CoreSpecimens and Determining Dimensional and Shape Tolerances,” ASTM D4435-01; “Standard Test Method for Elastic Moduli of Intact Rock CoreSpecimens in Uniaxial Compression,” ASTM D 3148-02; “Standard TestMethod for Triaxial Compressive Strength of Undrained Rock CoreSpecimens Without Pore Pressure Measurements,” ASTM D 2664-04; “StandardTest Method for Creep of Cylindrical Soft Rock Specimens in UniaxialCompressions,” ASTM D 4405-84, Reapproved 1989; “Standard Test Methodfor Performing Laboratory Direct Shear Strength Tests of Rock SpecimensUnder Constant Normal Stress,” ASTM D 5607-95; “Method of Test forDirect Shear Strength of Rock Core Specimen,” U.S. Military Rock TestingHandbook, RTH-203-80, available at“http://www.wes.army.mil/SL/MTC/handbook/RT/RTH/203-80.pdf” (lastaccessed on Jun. 25, 2010); and “Standard Method of Test for MultistageTriaxial Strength of Undrained Rock Core Specimens Without Pore PressureMeasurements,” U.S. Military Rock Testing Handbook, available athttp://www.wes.army.mil/SL/MTC/handbook/RT/RTH/204-80.pdf” (lastaccessed on Jun. 25, 2010). One of ordinary skill will recognize thatother methods of testing rock specimens may be used to determine thephysical constants used herein.

“Permeability” is the capacity of a rock to transmit fluids through theinterconnected pore spaces of the rock. Permeability may be measuredusing Darcy's Law: Q=(kΔP A)/(μL), where Q=flow rate (cm³/s),ΔP=pressure drop (atm) across a cylinder having a length L (cm) and across-sectional area A (cm²), μ=fluid viscosity (cp), and k=permeability(Darcy). The customary unit of measurement for permeability is themillidarcy. The term “relatively permeable” is defined, with respect toformations or portions thereof, as an average permeability of 10millidarcy or more (for example, 10 or 100 millidarcy). The term“relatively low permeability” is defined, with respect to formations orportions thereof, as an average permeability of less than about 10millidarcy. An impermeable layer generally has a permeability of lessthan about 0.1 millidarcy. By these definitions, shale may be consideredimpermeable, for example, ranging from about 0.1 millidarcy (100microdarcy) to as low as 0.00001 millidarcy (10 nanodarcy).

“Poisson's ratio” of an elastic material is a material property thatdescribes the amount of radial expansion when subject to an axialcompressive stress (or deformation measured in a direction perpendicularto the direction of loading). Poisson's ratio is the ratio of theelastic material radial deformation (strain) to its axial deformation(strain), when deformed within its elastic limit. Rocks usually have aPoisson's ratio ranging from 0.1 to 0.4. The maximum value of Poisson'sratio is 0.5 corresponding to an incompressible material (such aswater). It is denoted by the Greek letter ν (nu) and is a unit-lessratio.

A “reservoir” or “hydrocarbon reservoir” is defined as a pay zone (forexample, hydrocarbon-producing zones) that includes sandstone,limestone, chalk, coal, and some types of shale. Pay zones can vary inthickness from less than one foot (0.3048 m) to hundreds of feet(hundreds of m). The permeability of the reservoir formation providesthe potential for production.

“Reservoir properties” and “Reservoir property values” are defined asquantities representing physical attributes of rocks containingreservoir fluids. The term “Reservoir properties” as used in thisapplication includes both measurable and descriptive attributes.Examples of measurable reservoir property values include impedance top-waves, impedance to s-waves, porosity, permeability, water saturation,and fracture density. Examples of descriptive reservoir property valuesinclude facies, lithology (for example, sandstone or carbonate), andenvironment-of-deposition (EOD). Reservoir properties may be populatedinto a reservoir framework of computational cells to generate areservoir model.

A “rock physics model” relates petrophysical and production-relatedproperties of a rock formation to the bulk elastic properties of therock formation. Examples of petrophysical and production-relatedproperties may include, but are not limited to, porosity, pore geometry,pore connectivity volume of shale or clay, estimated overburden stressor related data, pore pressure, fluid type and content, clay content,mineralogy, temperature, and anisotropy and examples of bulk elasticproperties may include, but are not limited to, P-impedance andS-impedance. A rock physics model may provide values that may be used asa velocity model for a seismic survey.

“Shale” is a fine-grained clastic sedimentary rock with a mean grainsize of less than 0.0625 mm. Shale typically includes laminated andfissile siltstones and claystones. These materials may be formed fromclays, quartz, and other minerals that are found in fine-grained rocks.Non-limiting examples of shales include Barnett, Fayetteville, andWoodford in North America. Shale has low matrix permeability, so gasproduction in commercial quantities requires fractures to providepermeability. Shale gas reservoirs may be hydraulically fractured tocreate extensive artificial fracture networks around well bores.Horizontal drilling is often used with shale gas wells.

“Strain” is the fractional change in dimension or volume of thedeformation induced in the material by applying stress. For mostmaterials, strain is directly proportional to the stress, and dependsupon the flexibility of the material. This relationship between strainand stress is known as Hooke's law, and is presented by the formula:f/S=−Y (∂ξ/∂x).

“Stress” is the application of force to a material, such as a through ahydraulic fluid used to fracture a formation. Stress can be measured asforce per unit area. Thus, applying a longitudinal force f to across-sectional area S of a strength member yields a stress which isgiven by f/S.

The force f could be compressional, leading to longitudinallycompressing the strength member, or tensional, leading to longitudinallyextending the strength member. In the case of a strength member in aseismic section, the force will typically be tension.

“Thermal fractures” are fractures created in a formation caused byexpansion or contraction of a portion of the formation or fluids withinthe formation. The expansion or contraction may be caused by changingthe temperature of the formation or fluids within the formation. Thechange in temperature may change a pressure of fluids within theformation, resulting in the fracturing. Thermal fractures may propagateinto or form in neighboring regions significantly cooler than the heatedzone.

The “Young's modulus” of a rock sample is the stiffness of the rocksample, defined as the amount of axial load (or stress) sufficient tomake the rock sample undergo a unit amount of deformation (or strain) inthe direction of load application, when deformed within its elasticlimit. The higher the Young's modulus, the harder it is to deform, and,thus, the more brittle the rock. It is an elastic property of thematerial and is usually denoted by the English alphabet E having unitsthe same as that of stress.

Overview

Exemplary embodiments of the present techniques provide a systematicanalysis method to determine hydraulic fracture behavior of variousductile rocks. The method can employ an energy-based cohesive fracturemechanics model that incorporates unified creep-plasticity (UCP) into aconstitutive model. A temperature and strain-rate-dependent constitutiverelationship for ductile rock can be incorporated into the computationalmodel through finite element analysis. The singularity ahead of thefracture tip may be simulated using collapsed quadrilateral elements inthe finite element analysis. From the model that is developed, theimportant mechanical parameters for the simulation of ductile rockfracture behavior, such as the J-integral value and the stress fieldaround the fracture, can be numerically determined for hydraulicfracture modeling. The fundamental difference of linear elastic fracturemechanics and cohesive fracture mechanics-based analysis can then beexploited from both theoretical and computational aspects. Based on theanalysis of the fracture process zone of ductile rock, the effectivefracture material properties that account for the effect of rockductility are calculated and applied in currently available hydraulicfracture modeling tools.

FIG. 1 is a diagram 100 of a hydrocarbon reservoir 102, in accordancewith an exemplary embodiment of the present techniques. The hydrocarbonreservoir 102 may be a shale layer containing economically valuablequantities of natural gas. Wells 104, 106, and 108, may be drilled fromthe surface 110 through the overburden 112 to reach the reservoir 102.As indicated by the bidirectional arrows, each of the wells 104, 106,and 108 may be used for injection, for example, to fracture the rock ofthe reservoir 102, and for production, such as the harvesting ofhydrocarbons from the reservoir 102. The conversion from fracturing thereservoir 102 to harvesting hydrocarbons does not have to be limited toone cycle, as several cycles of fracturing and harvesting may bealternated.

The rock of the reservoir 102 may have a complex structure, for example,with faults 114 or other layers separating regions 116 and 118 of thereservoir 102. Accordingly, multiple wells may be used to access thedifferent regions. For example, in reservoir 102 shown in the diagram100, two wells 104 and 106 are used to access a first region 116 on oneside of a fault 114, while a third well 108 is used to access a secondregion 118 on the opposite side of the fault 114.

As mentioned previously, the relatively low permeability of rock in areservoir 102 may make fracturing the rock useful for harvestingeconomically viable amounts of hydrocarbons. The fracturing may beperformed by any number of techniques used for delivering energy to therock of the reservoir 102. For example, thermal fracturing may be usedto expand rock layers in or around the reservoir 102, causing stressesthat lead to the formation of fractures. Thermal fracturing may alsocause cracking by causing the expansion of fluids within the reservoir102, which creates hydraulic fractures. The present techniques are notlimited to any particular type of fracturing, as they are generallyfocused on the modeling of fracture growth in ductile rock layers.

Each of the wells 104, 106, or 108, may have perforations 120 through awell casing into the rocks of the reservoir 102, indicated by dots alongthe wells 104, 106, or 108. A fracturing fluid may be pumped at highpressure through the perforations into the reservoir 102. When thepressure of the fluid exceeds the strength of the rock layer, fractureswill form as the rock breaks. The fracturing fluid may be aqueous ornon-aqueous, and may contain materials to hold open fractures that form.These materials, known as proppants, may include natural or syntheticmaterials such as sand, gravel, ground shells, glass beads, or metalbeads, among others. The hydraulic fracturing process is discussed infurther detail with respect to FIG. 2.

FIG. 2 is a diagram 200 that shows hydraulic fracturing through a singlewell 202 extending from a surface facility 204 through an overburden 206into a reservoir 208, in accordance with an exemplary embodiment of thepresent techniques. Perforations along the well 202 can allowpressurized fracturing fluid to flow into the formation, creating afracture zone 210 around the perforated section of the well 208. Thefracture zone 210 may be considered a network or “cloud” of fracturesgenerally radiating out from the well 202. Each fracture is a crackformed in the rock, as discussed with respect to a close up view shownin FIG. 3.

FIG. 3 is a close-up view of a fracture 302 forming in the reservoir 208of FIG. 2. As shown in FIG. 3, a zone 304 may be created ahead of thetip of the fracture 302. Previous models, based on linear-elasticfracture mechanics (LEFM) often assumed that the zone 304 at the tip ofthe fracture 302 had minimal volume. This assumption provided inaccurateresults for some materials, as discussed further with respect to FIG. 4.

FIG. 4 shows a non-linear fracture zone 304 ahead of a fracture 302 inmaterials ranging from brittle to ductile. The non-linear fracture zone304 ahead of the fracture 302 consists of two zones: the fractureprocess zone 402 characterized by progressive softening and a nonlinearhardening zone 404 characterized by perfect plasticity or plastichardening. The stress versus distance in front of each crack tip isshown in the graphs on the right. FIG. 4( a) shows a non-linear fracturezone 304 for a very brittle material, such as glass, ceramics, and thelike. Since the entire non-linear fracture zone 304 is small (e.g., asingularity) compared to the size of the fracture 302, linear elasticfracture mechanics (LEFM) is applicable in this case, and the materialtoughness term, K, dominates. FIG. 4( b) shows the non-linear fracturezone 304 for a very ductile material where the softening or fractureprocess zone 402 is still small, but the nonlinear hardening zone 404 isnot negligible, such as ductile metals, tough alloys, many plastics, andthe like. In this case, the nonlinear hardening zone 404 dominates thefracture behavior and elasto-plastic fracture mechanics using the energyof the fracture, e.g., the J-integral, can be employed.

FIG. 4( c) shows the behavior of a quasi-brittle material, such asconcrete, cemented aggregates, many rocks, and the like. In thesematerials, a major part of the non-linear fracture zone 304 is afracture process zone 402, which undergoes progressive damage withmaterial softening due to micro-cracking and void formation. Thecracking in the fracture process zone can be distributed over a finitesize volume, which exhibits so-called strain-softening, for example, astress strain relation in which the maximum principal stress decreaseswith corresponding increasing principal strain. Here, cohesive fracturemechanics-based theory can be applied, in accordance with exemplaryembodiments of the present techniques. The ratio of the fracture processzone 402 to the nonlinear zone 304 in a quasi-brittle material may beabout 1:1, 1:1.2, or 1:1.5.

Most ductile rock, such as ductile shale, follows this behavior. Inexemplary embodiments of the present techniques, a constitutive modelfor rock behavior is revised to incorporate plasticity and creep effectsin modeling of hydraulic fractures, as discussed with respect to FIG. 5,creating a UCP constitutive model. A cohesive zone model or CZM can thenbe used to determine the fracture properties.

Incorporation of Unified Creep Plasticity Theory into a ConstitutiveModel to Simulate Temperature and Strain-Rate-Dependent Behavior ofDuctile Rock.

FIG. 5 is a process flow chart of a method 500 for modeling fractures inductile rock, in accordance with an exemplary embodiment of the presenttechniques. The method begins at block 502 with the generation of aconstitutive model that incorporates a unified creep-plasticity (UCP)theory.

In an exemplary embodiment of the present techniques, an appropriaterock constitutive model interrelates the behavior of ductile rock, suchas shale, with load and boundary conditions. As noted above, the modeldiffers from models for more brittle rocks by incorporating creep andplasticity effects into the constitutive model for ductile rock. Both ofthese effects, creep and plasticity, arise in ductile rock during longterm, high rate injection or production of fluids. Unifiedcreep-plasticity (UCP) theory can be used to address these physicaleffects.

The constitutive model may be strain rate and temperature dependent. TheUCP theory can be incorporated into the finite element model based onthe viscoplasticity theory, as would be understood by one of skill inthe art. The model can be programmed into a finite element modelingsoftware. For example, a user-defined subroutine may be linked to thefinite element analysis (FEA) analysis package. The theory givesaccurate predictions for the constitutive relationship of ductile rock.Generally, some experimentally determined material constants are neededto describe the mechanical behaviors and damage evolution to achieve areasonable accuracy, such as tensile strength, Young's modulus,Poisson's ratio, fracture toughness etc. The stress-strain relationshipcan be measured from a standard tri-axial or poly-axial experiment, andother material parameters, such as the A and B discussed with respect toEqn. 2, can be fitted using experimental data. Such as tensile strength,Young's modulus, Poisson's ratio, fracture toughness etc. Thestress-strain relationship can be measured from a standard tri-axial orpoly-axial experiment, the other material parameters such as A, B can befitted using experimental data.

For low or impermeable ductile rock, the relationship between stressrate a and strain rate can be written as shown in Eqn. 1.

$\begin{matrix}{\overset{.}{\sigma} = {{C(T)}:{( {\overset{.}{ɛ} - {\overset{.}{ɛ}}^{in} - {\phi \; T\; I}} ) + {C^{- 1}:{\sigma:{\frac{\partial C}{\partial T}\overset{.}{T}}}}}}} & {{Eqn}.\mspace{11mu} 1}\end{matrix}$

In Eqn. 1, C represents the matrix of the elastic constants and is afunction of Poisson's ratio ν, Young's modulus E(T), and absolutetemperature T. The constants represents the total strain rate; {dot over(ε)}^(in) represents the inelastic strain rate; φ represents thecoefficient of thermal expansion and I represents the unit tensor.Accordingly, the term φT I refers to the thermal strain rate.

The form of the flow rule for inelastic strain rate can be given by theformula in Eqn. (2):

$\begin{matrix}{{\overset{.}{ɛ}}^{in} = {\sqrt{\frac{3}{2}}{A( \frac{S_{v}}{d} )}^{n}{\exp ( {B( \frac{S_{v}}{d} )}^{n + 1} )}{\Theta (T)}N}} & {{Eqn}.\mspace{11mu} 2}\end{matrix}$

In Eqn. 2, d represents the drag strength, A and B represent materialconstants; and the diffusivity term Θ(T)=exp(−Q/(kT)). Here A is a creepmaterial constant, B is just a unitless material constant, both of theparameters can be fit from experimental data. Q represents the apparentactivation energy and k is the universal gas constant.N represents a unit vector in the direction of the deformation loadingdefined by: N=(S−α)/∥S−α∥, where the deviatoric stress tensor can beobtained by subtracting the hydrostatic stress tensor from the stresstensor defined by: S=σ−tr(σ)/3; S_(ν) represents the viscous overstressdefined as

${S_{v} = {\langle{{\sqrt{\frac{3}{2}}{{S - \alpha}}} - R}\rangle}},$

in which R is the yield stress radius and α is the deviatoric backstress. The evolutionary law for α is defined as {dot over (α)}μ∥{dotover (ε)}^(in)∥N−βα, where μ and β are material constants. It can beunderstood that α reaches a saturated value along with the loadingprocess. The R term is determined by R=R_(initial)+R_(add) and {dot over(R)}_(add)= ω{dot over (χ)}, in which ω is a material constant and χ isthe hardening from the dislocation increment within the material, whichis defined as {dot over (χ)}=μ∥{dot over (ε)}^(in)∥( χ−χ), where χ isthe sutured value of χ; μ is a material constant. The material constantsfor the model can be obtained through analysis of the material behaviorand by fitting to the experimental data for different rocks. Suchmaterial constants may include tensile strength, fracture toughness,Young's modulus, permeability, A, and B, among others.

The result of block 502 is a UCP constitutive model that may be used topredict fractures in quasi-brittle materials, for example, as discussedwith respect in FIG. 4( c). Such materials may include shales,concretes, plastic composites, and other materials having a largefracture process zone.

Predicting the J-Integral Value and the Stress Field from Finite ElementAnalysis, Using a UCP Constitutive Model that Accounts for Ductility.

At block 504, the constitutive model for ductile rock defined at block502 can be analyzed by performing finite element analysis to predict thevalue of the J-integral and the stress field considering plasticdeformation energy. To simulate the singularity at the fracture tip, aring of collapsed quadrilateral elements can be used to perform themodeling in the fracture tip region, as discussed further with respectto FIG. 6.

FIG. 6 is a graph 600 of a von Mises stress distribution ahead of afracture tip, in accordance with an exemplary embodiment of the presenttechniques. To simulate the singularity at the tip 602 of the fracture604, a ring 606 of collapsed quadrilateral elements can be used toperform the modeling in a region 608 surrounding the tip 602 of thefracture 604. For example, in an exemplary embodiment, the quadrilateralelements may be the CPS8R elements in the finite element analysisprogram ABAQUS, available from Dassault Systèmes under the SIMULIAbrand. The mid-side nodes 610 on the sides, connected to the fracturetip, can be moved to the ¼ point nearest to the fracture tip, and aquarter point spacing with second-order isoparametric elements can thenbe created. The J-integral and stress field can be predicted numericallyfrom the model by considering plasticity effect combined with revisedrock constitutive model.

More specifically, the J-integral value and the stress field can bepredicted from the finite element analysis by employing the UCPconstitutive model to consider ductility, e.g., the plastic deformationenergy. The three-dimensional expression of the J-integral can bedefined as shown in Eqn. 3.

$\begin{matrix}{J_{l} = {{{- {\int_{\Gamma}{\sigma_{ij}n_{j}u_{i,l}\ {S}}}} + {\int_{\Gamma}{W\; n_{l}\ {S}}}} = {\int_{\Gamma}{P_{lj}n_{j}\ {S}}}}} & {{Eqn}.\mspace{14mu} 3}\end{matrix}$

In Eqn. 3, the Maxwell energy-momentum tensor isP_(lj)=Wδ_(lj)−σ_(ij)u_(i,l) and W=∫₀ ^(ε) ^(ij) σ_(ij)dε_(ij) andW_(,l)=σ_(ij)ε_(ij,l). Further,

$ɛ_{ij} = {\frac{1}{2}( {u_{i,j} + u_{j,i}} )}$

and the integration path Γ is an arbitrary closed contour followedcounterclockwise surrounding the fracture tip in a stressed body.Additional, in Eqn. 3, σ_(ij) represents the stress, ε_(ij) representsthe strain, u_(i) represents the displacement, and n_(j) represents thenormal vector on the boundary dS.

It may be assumed that the released energy rate to create a crack orfracture is equal to the value of the J-integral,

${J = {G = \frac{K_{c}^{2}}{E^{\prime}}}},$

in which J represents the J-integral value, G represents the energyrelease rate, K_(c) represents the fracture toughness and E′ representsthe Young's modulus. The fracture toughness is equal to √{square rootover (JE′)}, and the value of J is determined from numerical simulation,which accounts for plastic deformation and energy as discussed withrespect to block 502.

Because the plasticity deformation and energy has been incorporated inthe finite element analysis through the revised rock constitutive model,the effective fracture toughness value can be determined numericallyfrom the predicted J-integral value. The effective fracture parametersmethod can be applied to consider plasticity effect in hydraulicfracture, as discussed below.

Determination of the Ductile Rock Fracture Parameters by CohesiveFracture Mechanics and Numerical Modeling.

Returning to FIG. 5, the effective fracture properties for the ductilerock may be determined at block 506 by either theoretical or numericalmethods. One of the reasons that LEFM-based methods may fail to giveaccurate predictions for ductile rock is that these methods neglect thefracture process zone ahead of the fracture tip, which may besignificant in ductile rock fracture analysis. In LEFM, R is oftenassumed to be infinitesimally small.

Theoretically, the fundamental mechanical difference of LEFM andcohesive fracture mechanics-based methods can be explained by the sizeof the fracture process zone and its effect on fracture extension inductile rock. The softening behavior ahead of the fracture tip can beconsidered by employing revised fracture properties. For ductile shale,which is quasi-brittle in nature, the effective fracture toughness canbe derived to accommodate the ductile behavior based on the fracturecharacteristic length ahead of fracture tip, as discussed with respectto FIG. 7.

FIG. 7 is graph illustrating the estimation of the fracture toughnessfor a nonlinear softening material-ductile rock based on the size of thefracture process zone, in accordance with an exemplary embodiment of thepresent techniques. Based on size of the fracture process zone, R_(C),ahead of tip 702 of the fracture 704, a straightforward method can bedeveloped to consider the effect of ductility on fracture parameters. Inthis graph 700, the y-axis 706 represents stress, while the x-axis 708represents the distance from the fracture tip 702. The highest stresspoint 710 has a stress represented by σ_(y), and occurs at the tip ofthe fracture process zone 304 that is farthest from the fracture tip702, as discussed with respect to FIG. 4. By assuming a finite size,fracture process zone, R_(C), ahead of the fracture tip 702, the lengthof the nonlinear fracture zone and softening of the ductile rock can bedetermined and used to estimate the fracture toughness. The rocksoftening may be assumed to be parabolic, resulting in a formula for thestress of a first curve 712, located before the end of the fractureprocess zone R_(C) of

$\sigma_{yy} = {{f_{t}( \frac{x}{R_{C}} )}^{n}.}$

A second curve 714 after the fracture process zone R_(C) may have aformula for the stress of

$\sigma_{yy} = {\frac{K_{IC}}{\sqrt{2{\pi ( {x - r_{1}} )}}}.}$

In these equations, n represents the parabolic distribution degree,σ_(yy) represents the calculated stress, f_(t) represents rock tensilestrength, and x represents the distance from the fracture tip 702. R_(C)represents the length of the fracture process zone, r₁ represents thedistance of the equivalent elastic crack to be shifted ahead of the truefracture tip, and K_(IC) represents the effective fracture toughness ofthe rock.

The y-axis 706 may be shifted by drawing a new y-axis, y′ 716 at thepoint at which the area under the second curve 714 from x=r₁ to R_(C) isequal to the area under the first curve 712 from x=0 to R_(C). Applyingthis “equivalent area rule” to the equations above provides the formulashown in Eqn. 4.

$\begin{matrix}{{\int_{r\; 1}^{Rc}{{K_{IC}\lbrack {2{\pi ( {x - r_{1}} )}} \rbrack}^{{- 1}/2}\ {x}}} = {{\int_{0}^{Rc}{{f_{t}^{\prime}( \frac{x}{R_{c}} )}^{n}{x}}} = {\frac{1}{n + 1}R_{c}f_{t}^{\prime}}}} & {{Eqn}.\mspace{14mu} 4}\end{matrix}$

After integration, Eqn. 4 gives

${R_{c} = {{\frac{n + 1}{\pi}( \frac{K_{IC}}{f_{t}^{\prime}} )^{2}} = ( \frac{K_{eff}}{f_{t}^{\prime}} )^{2}}},$

and the effective fracture toughness can be defined by the formula shownin Eqn. 5.

$\begin{matrix}{K_{eff} = {\sqrt{\frac{n + 1}{\pi}}K_{IC}}} & {{Eqn}.\mspace{14mu} 5}\end{matrix}$

In these formulas, the value of the softening parabolic distributiondegree n can be determined for different rock materials to include thefracture tip effect on fracture growth. The idealized characteristiclength of fracture can be defined by

$l_{c} = {( \frac{K_{IC}}{f_{t}^{\prime}} )^{2}.}$

The actual characteristic length of the fracture can then be representedby

$l_{c} = {\frac{n + 1}{\pi}{( \frac{K_{IC}}{f_{t}^{\prime}} )^{2}.}}$

For typical quasi-brittle materials such as concrete, the value of

$\frac{n + 1}{\pi}$

usually varies from 2 to 5, for example, K_(eff)=1.414˜2.236_(IC). Forductile rock, the value of n can be experimentally or numericallydetermined from the value of the J-integral considering plastic straineffect to the fracture energy by fitting experimental or numericalanalysis data to a rock stress-strain relationship. Thus, the numericalmodel proposed herein can be applied to determine the effective fractureparameter because plasticity strain is incorporated in the rockconstitutive model. The numerical investigation of effective materialproperties such as fracture toughness can be determined from the resultsof computational model discussed herein.

Predict Hydraulic Fracturing of Ductile Rock Based on the Results fromPrevious Steps.

Returning to FIG. 5, at block 508, the model developed above can be runto predict fracture formation in ductile rock, for example, fromhydraulic fracturing. In an exemplary embodiment, a pore pressurecohesive zone model (pore pressure CZM) is developed to predicthydraulic fracture during injection using the model developed above. Thepore pressure CZM is a numerical approach developed to model fractureinitiation and growth that may be used in quasi-brittle materials thathave a material softening effect. The method 500 treats the fracture asa gradual process in which separation between incipient materialsurfaces is resisted by cohesive tractions. The fracture initiation andpropagation can be predicted from the pore pressure CZM analysis throughthe application of a finite element analysis incorporating the UCPconstitutive model.

The damage mechanics incorporated in the model can be introduced and thedamage initiation and evolution conditions can be defined. The porepressure CZM has been applied herein to hydraulic fracture growth usingdifferent rock properties. Based on the computational and theoreticalanalysis discussed above, a revision to fracture parameters, such asfracture toughness, can be implemented to account for the ductileeffects.

In an exemplary embodiment, the pore pressure CZM is used to predicthydraulic fracture of injection wells. By combining CZM with the UCPconstitutive model, discussed with respect to block 502, the hydraulicfracture of ductile rock can be predicted directly though finite elementanalysis. For other hydraulic fracture model or CZM with an elasticconstitutive model, the effective fracture parameter method can beapplied to predict plasticity effect conveniently.

To incorporate a cohesive zone model into finite element analysis, theprinciple of virtual work can be modified to the formula shown in Eqn.6.

∫_(V) {tilde over (P)}:δ{tilde over (F)}dV−∫ _(S) _(i) {right arrow over(T)} _(CZ) ·δΔ{right arrow over (u)}dS=∫ _(δV) {right arrow over (T)}_(e) ·δ{right arrow over (u)}dS  Eqn. 6

The principle assumes that the total work done by all forces acting on asystem in static equilibrium is zero for any infinitesimal displacementfrom equilibrium, which is consistent with the constraints of thesystem. This may also be known as the virtual work principle. In Eqn. 6,{tilde over (P)} represents the nominal stress tensor, which equals{tilde over (F)}⁻¹ det({tilde over (F)}){tilde over (σ)}. {tilde over(F)} represents the deformation gradient and {tilde over (σ)} representsthe Cauchy stress tensor. V represents the control volume, and {rightarrow over (u)} represents the displacement vector. {right arrow over(T)}_(CZ) represents the cohesive zone traction vector and {right arrowover (T)}_(e) represents the traction vector on the external surface ofthe body. S_(i) represents the internal boundary and ∂V represents theexternal boundary of V.

FIG. 8 is a graph 800 illustrating the cohesive elementtraction-separation (T-S) Law used in the CZM, in accordance with anexemplary embodiment of the present technique. The Traction-separation(T-S) Law may be used to define the constitutive response of a cohesiveelement layer directly in terms of traction (applied force) versusseparation (fracture opening/displacement). In FIG. 8, the y-axis 802represents the tensile force (T) and the x-axis 804 represents theseparation distance (δ). The maximum height is the maximum tensileforce, T_(ult), reached before a fracture initiates at δ_(ini). Thepoint at which no tensile forces remain on the system is δ_(fail), whichrepresents complete failure of the rock. The area 806 under the curverepresents the energy placed on the rock. The area 808 under the curveup to δ_(ini), represents the energy needed to initiate failure.

The leak-off model in pore pressure CZM is based on hydraulic fluidcontinuity. The tangential permeability, k_(t) is defined based onReynold's equation, k_(t)=d³/12μ, where μ represents the fluidviscosity, d represents the gap opening defined byd=t_(curr)−t_(orig)+g_(init), and t_(curr) and t_(orig) are the currentand original cohesive element geometrical thicknesses, respectively. Aswould be understood, in the ABAQUS software, it is actually resistanceto fluid flow rather than permeability. The term g_(init) is the initialgap opening, which has a default value of 0.002. Permeability ofreservoir rock, or hydraulic conductivity, k, is used in a flow modelbased on Darcy's equations for permeability

$\hat{k} = {{\frac{v}{g}\overset{\_}{k}} = {\frac{\mu}{\rho \; g}\overset{\_}{k}\mspace{14mu} {( {m\text{/}s} ).}}}$

In this equation, ν represents kinematic viscosity, μ represents dynamicviscosity, ρ represents fluid density (=1 for water), and g representsgravity acceleration (=9.8 m/s²).

A quadratic nominal stress criterion can be used to predict damageinitiation in the pore pressure CZM. Damage can be assumed to initiatewhen a quadratic interaction function involving the nominal stressratios reaches unity. This criterion can be represented as

${( \frac{\langle t_{n}\rangle}{t_{n}^{0}} )^{2} + ( \frac{t_{s}}{t_{s}^{0}} )^{2} + ( \frac{t_{t}}{t_{t}^{0}} )^{2}} = 1.$

In this equation, t_(n), t_(s), t_(t) refer to the normal, the first,and the second shear stress components. The terms, t_(n) ⁰, t_(s) ⁰,t_(t) ⁰ represent the peak values of the nominal stress when thedeformation is either purely normal to the interface or in the first orthe second shear direction. The symbol

is the usual Macaulay bracket interpretation, used to signify that apurely compressive deformation or stress state does not initiate damage.

The damage evolution for mixed mode failure is defined based onBenzeggagh-Kenane fracture criterion, when the critical fractureenergies during deformation along the first and the second sheardirections are similar; i.e. G_(s) ^(C)=G_(t) ^(C). The criteria can begiven by

${{G_{n}^{C} + {( {G_{s}^{C} - G_{n}^{C}} )( \frac{G_{S}}{G_{T}} )^{\eta}}} = G^{C}},$

in which the mixed-mode fracture energy, G^(C)=G_(n)+G_(s)+G_(t). G_(n),G_(s), and G_(t) refer to the work done by the traction and itsconjugate relative displacement in the normal, the first, and the secondshear directions, respectively. G_(n) ^(C), G_(s) ^(C) and G_(t) ^(C)refer to the critical fracture energies required to cause failure in thenormal, the first, and the second shear directions, respectively. η is amaterial parameters that represents softening behavior, i.e., a shape ofthe softening curve, and may be determined by polyaxial or triaxial rockcell experiments.

The effective material properties predicted from the use of the porepressure CZM can also be used in other current available LEFM-basedhydraulic fracture tools. For example, the effective fracture toughnessmay be used in the pseudo 3D fracture model to revise the prediction byconsidering ductility effects.

Examples

A short-term water injection case with zero leak-off was analyzed toinvestigate the accuracy of the fracture geometry prediction using apore pressure CZM based on the UCP constitutive model described herein.The predicted fracture length and δ_(fail) is given in Table 1. Theresult shows that, compared with previous LEFM-based models, forexample, a standard Pseudo 3D model and the Perkins-Kern-Nordgrenmethod, the pore pressure CZM-based method predicted the fracture lengthmore accurately.

TABLE 1 Fracture length predicted using different models AnalyticalPseudo solution (Dean CZM 3D PKN et al. 2008) Fracture 147.5 180 183141.5 length (ft) δ_(fail) (in) 0.0059 — — 0.0060

As another example, the pore pressure CZM was used to predict theeffects of different parameters on fracture geometry, as shown in Table2. Based on the parametric analysis, increases in Young's modulus (E) ofthe rock results in decrease of the fracture width, while the fracturelength and height increases. Increase of injection rate, increasesfracture width, length and height. As shown in Table 2 and discussedwith respect to FIG. 13, increasing the leak-off coefficient decreasesfracture width, length, and height. Increasing injected fluid viscosityresults in decreasing of the fracture length, while fracture width andheight increase.

TABLE 2 Effects of different parameters to fracture geometry predictedusing pore pressure CZM Fracture Fracture Fracture Parameter widthLength Height Rock E ↑ ↓ ↑ ↑ Inj. rate ↑ ↑ ↑ ↑ Leak-off coef. ↑ ↓ ↓ ↓Fracture fluid viscosity ↑ ↑ ↓ ↑

FIG. 9 is a drawing of a three-layer 3D pore pressure cohesive elementmodel 900 that may be used to predict hydraulic fracture caused by fluidinjection using exemplary embodiments of the current techniques. Theoverburden 902, side burden 904, and pore pressure 906 boundaryconditions are applied. The model 900 has three layers, a top layer 908of hard rock, a middle layer 910, and a lower layer 912 of ductile rock.The top and bottom layers 908 and 912 have perforations (or fractures)914, for example, caused by pressure from the injection of a fracturingfluid, as indicated by larger arrows 916. The results from this modelare discussed with respect to FIG. 10.

FIG. 10 is an illustration of a fracture geometry predicted for themodel of FIG. 9, in accordance with exemplary embodiments of the presenttechniques. A pore pressure CZM was used to predict a first fracture1002 into the brittle top layer 908 using an elastic (LEFM) constitutivemodel. The pore pressure CZM was used to predict a second fracture 1004into the quasi-brittle bottom layer 912 using a plastic constitutivemodel, as described herein. It can be seen that the first fracture 1002,in the hard rock of the top layer 908, is longer and narrower than thesecond fracture 1004 in the ductile rock of the bottom layer 912. Ashorter fracture 1004 in the ductile rock layer 912 may be because alarger deformation in a ductile rock absorbs more energy than in a hardrock, and, therefore, less energy is being used to extend the secondfracture 1004.

In another example, for a pseudo 3D hydraulic model or finite elementanalysis with elastic rock constitutive model, the effective fractureparameter can be used to simulate the plasticity effect by applying themethod developed in the present techniques. FIGS. 11 and 12 illustrate afracture propagation prediction with different effective fracturetoughness values using the finite element analysis and pseudo 3D model,respectively.

FIG. 11 is a series of graphs illustrating predicted fracturepropagation prediction of finite element analysis using a UCPconstitutive model, in accordance with exemplary embodiments of thepresent techniques. The finite element analysis shown in each of thefour graphs was performed at different calculated fracture toughnessvalues. Thus, to generate the graph of FIG. 11( a), the fracturetoughness, K_(eff), was set equal to the toughness under the LEFMassumptions, K_(Ic). To generate FIG. 11( b), K_(eff) was set to 1.5K_(Ic). To generate FIG. 11( c), K_(eff) was set to 2 K_(Ic). Finally,to generate FIG. 11( d), K_(eff) was set to 3 K_(Ic). As shown in thesegraphs, the fracture length decreases as the effective fracturetoughness increases, simulating the plasticity effect. This matches thephysical observation made in field applications. Further confirmation ofthis effect can be provided using other models, for example, asdiscussed with respect to FIG. 12.

FIG. 12 is a series of graphs from a finite element analysis thatillustrate changes in the fracture propagation predictions from a pseudo3D model, in accordance with exemplary embodiments of the presenttechniques. Each of the graphs uses the same K_(eff) ratio to K_(Ic) asthe corresponding graph in FIG. 11. Thus, to generate FIG. 12( a)K_(eff) was set equal to K_(Ic). To generate FIG. 12( b), K_(eff) wasset to 1.5 K_(Ic). To generate FIG. 12( c), K_(eff) was set to 2 K_(Ic).Finally, to generate FIG. 12( d), K_(eff) was set to 3 K_(Ic). Theresults confirm those shown in FIG. 11. As the effective fracturetoughness increases the effective fracture length decreases.

Further tests were run using the model to determine the effects ofincreasing the leak-off coefficient on the fractures. As notedpreviously, a higher leak-off coefficient decreases the effect offracture toughness on normalized fracture length. For higher leak-offcoefficient formation, the leak-off phenomenon dominates the hydraulicfracture process, while the fracture toughness plays a more importantrole with a higher fracture efficiency of injected fluid. Fractureefficiency is defined as the ratio of the fracture volume to the volumeof fluid injected. For low viscosity fluids, fracture toughness can beone of the dominant parameters controlling fracture growth as discussedwith respect to FIG. 13.

FIG. 13 is a graph 1300 of effective fracture toughness at differentleak-off coefficients, in accordance with exemplary embodiments of thepresent techniques. From the graph 1300, it can be seen that for a highpermeability formation, for example, with leak-off coefficients >0.005,the fracture toughness has less effect on fracture length. Therefore,the application of effective fracture toughness is not necessary.However, for low permeability ductile rock with leak-off coefficientsless than 0.001, effective fracture toughness has significant effect onfracture geometry. In embodiments, effective fracture toughness may beused to simulate ductile effect in hydraulic fracture.

Conclusion

In exemplary embodiments, a combined theoretical and numerical method,based on cohesive fracture mechanics and unified creep and plasticitytheory, has been used to predict the plasticity effect in hydraulicfracture. An effective fracture parameter for quasi-brittle/ductile rockis proposed. Pore pressure cohesive zone model has been applied toinvestigate the hydraulic fracture behavior. With the revised UCPconstitutive model and the energy based J-integral calculation, therevised fracture toughness can be calculated to include the ductilityeffect into rock hydraulic fracture. The example indicates that for lowpermeability ductile rock, the present techniques can predict hydraulicfracture behavior more accurately than currently available modelingtools for hydraulic fracture.

Computing System

FIG. 14 is a block diagram of an exemplary cluster computing system 1400that may be used to implement the present techniques, in accordance withexemplary embodiments of the present techniques. The cluster computingsystem 1400 illustrated has four computing units 1402, each of which mayperform calculations for a part of the finite element analysiscalculations. However, one of ordinary skill in the art will recognizethat the present techniques are not limited to this configuration, asany number of computing configurations may be selected. For example, asmaller analysis may be run on a single computing unit 1402, such as aworkstation, while a large finite element analysis calculation may berun on a cluster computing system 1400 having 10, 100, 1000, or evenmore computing units 1402.

The cluster computing system 1400 may be accessed from one or moreclient systems 1404 over a network 1406, for example, through a highspeed network interface 1408. The computing units 1402 may also functionas client systems, providing both local computing support, as well asaccess to the wider cluster computing system 1400.

The network 1406 may include a local area network (LAN), a wide areanetwork (WAN), the Internet, or any combinations thereof. Each of theclient systems 1404 may have non-transitory, computer readable memory1410 for the storage of operating code and programs, including randomaccess memory (RAM) and read only memory (ROM). The operating code andprograms may include the code used to implement all or any portions ofthe methods discussed herein, for example, as discussed with respect toFIG. 5. Further, the non-transitory computer readable media may holdconstitutive models, unified creep theory models, and a combined UCPconstitutive model as discussed above. The client systems 1404 can alsohave other non-transitory, computer readable media, such as storagesystems 1412. The storage systems 1412 may include one or more harddrives, one or more optical drives, one or more flash drives, anycombinations of these units, or any other suitable storage device. Thestorage systems 1412 may be used for the storage of UCP constitutivemodels, code, data, and other information used for implementing themethods described herein.

The high speed network interface 1408 may be coupled to one or morecommunications busses in the cluster computing system 1400, such as acommunications bus 1414. The communication bus 1414 may be used tocommunicate instructions and data from the high speed network interface1408 to a cluster storage system 1416 and to each of the computing units1402 in the cluster computing system 1400. The communications bus 1414may also be used for communications among computing units 1402 and thestorage array 1416. In addition to the communications bus 1414 a highspeed bus 1418 can be present to increase the communications ratebetween the computing units 1402 and/or the cluster storage system 1416.

The cluster storage system 1416 can have one or more non-transitory,computer readable media devices, such as storage arrays 1420 for thestorage of UCP constitutive models, data, visual representations,results, code, or other information, for example, concerning theimplementation of and results from the methods of FIG. 5, such as, forexample, the finite element analysis results. The storage arrays 1420may include any combinations of hard drives, optical drives, flashdrives, holographic storage arrays, or any other suitable devices.

Each of the computing units 1402 can have a processor 1422 andassociated local tangible, computer readable media, such as memory 1424and storage 1426. Each of the processors 1422 may be a multiple coreunit, such as a multiple core CPU or a GPU. The memory 1424 may includeROM and/or RAM used to store code, for example, used to direct theprocessor 1422 to implement the methods discussed with respect to FIG.5. The storage 1426 may include one or more hard drives, one or moreoptical drives, one or more flash drives, or any combinations thereof.The storage 1426 may be used to provide storage for UCP constitutivemodels, intermediate results, data, images, or code associated withoperations, including code used to implement the method of FIG. 5.

In other aspects the present techniques and claimed subject matter mayinclude for foreign filing purposes:

1. A method (500) for modeling fractures (302), comprising generating(502) a model that incorporates a unified creep-plasticity (UCP)representation into a constitutive model for a ductile rock (208).

2. The method (500) of paragraph 1, wherein the ductile rock (208) formsfractures (302) that have a fracture process zone (402), wherein thefracture process zone (402) is distributed over a finite size volume andwherein the fracture process zone (402) exhibits strain-softening.

3. The method (500) of paragraph 2, wherein a ratio of a length of thefracture process zone (402) to a length of a non-linear fracture zone(304) located at the tip of a fracture (302) is between about 1:1 and1:1.5.

4. The method (500) of paragraph 1, further comprising predicting (504)a J-integral value and a stress field that incorporates a plasticityeffect from finite element analysis using the model.

5. The method (500) of paragraph 1, further comprising determining (506)an effective fracture property of the ductile rock (208).

6. The method (500) of paragraph 5, further comprising determining theeffective fracture property of the ductile rock (208) using a numericalmethod.

7. The method (500) of paragraph 1, further comprising predicting (508)a fracture (302) in the ductile rock (208).

8. The method (500) of paragraph 7, further comprising predicting ahydraulic fracture (302) of the ductile rock (208) using a pore pressurecohesive zone model (CZM).

9. The method (500) of paragraph 7, further comprising determining aflow of hydrocarbon in a fractured ductile rock (208).

10. The method (500) of paragraph 7, further comprising adjusting aparameter in a fracturing process based, at least in part, on thepredicted fracture (302) in the ductile rock (208).

11. The method of paragraph 10, wherein the parameter is a mass flowrate of a fluid, a pressure of a fluid, or a combination thereof.

12. A system (1400) for modeling fractures (302) in a reservoir (102),comprising: a processor (1402); a storage medium (1412, 1420, 1426)comprising a data structure comprising a model that incorporates unifiedcreep-plasticity (UCP) into a constitutive model for a ductile rock; anda computer readable medium (1410, 1424, 1412, 1420, 1426) comprisingcode configured to direct the processor (1402) to perform a finiteelement analysis of the ductile rock (208) using the model.

13. The system of paragraph 12, wherein the processor (1402) comprises amulti-processor cluster computer system (1400).

14. The system of paragraph 12, wherein the computer readable medium(1410, 1424, 1412, 1420, 1426) comprises code configured to direct theprocessor (1402) to run a finite element analysis for predicting aJ-integral value and a stress field incorporating plasticity effectusing the model.

15. The system of paragraph 12, wherein the computer readable mediumcomprises (1410, 1424, 1412, 1420, 1426) code configured to direct theprocessor (1402) to predict a hydraulic fracture (302) in the ductilerock (208).

The present techniques are not limited to the architecture or unitconfiguration illustrated in FIG. 14. For example, any suitableprocessor-based device may be utilized for implementing all or a portionof embodiments of the present techniques, including without limitationpersonal computers, networks personal computers, laptop computers,computer workstations, GPUs, mobile devices, and multi-processor serversor workstations with (or without) shared memory. Moreover, embodimentsmay be implemented on application specific integrated circuits (ASICs)or very large scale integrated (VLSI) circuits. In fact, persons ofordinary skill in the art may utilize any number of suitable structurescapable of executing logical operations according to the embodiments.

While the present techniques may be susceptible to various modificationsand alternative forms, the exemplary embodiments discussed above havebeen shown only by way of example. However, it should again beunderstood that the present techniques are not intended to be limited tothe particular embodiments disclosed herein. Indeed, the presenttechniques include all alternatives, modifications, and equivalentsfalling within the true spirit and scope of the appended claims.

What is claimed is:
 1. A method for modeling fractures, comprisinggenerating a model that incorporates a unified creep-plasticity (UCP)representation into a constitutive model for a ductile rock.
 2. Themethod of claim 1, wherein the ductile rock forms fractures that have afracture process zone, wherein the fracture process zone is distributedover a finite size volume, and wherein the fracture process zoneexhibits strain-softening.
 3. The method of claim 2, wherein a ratio ofa length of the fracture process zone to a length of a non-linearfracture zone located at the tip of a fracture is between about 1:1 and1:1.5.
 4. The method of claim 1, further comprising predicting aJ-integral value and a stress field that incorporates a plasticityeffect from finite element analysis using the model.
 5. The method ofclaim 1, further comprising determining an effective fracture propertyof the ductile rock.
 6. The method of claim 5, further comprisingdetermining the effective fracture property of the ductile rock using anumerical method.
 7. The method of claim 1, further comprisingpredicting a fracture in the ductile rock.
 8. The method of claim 7,further comprising predicting a hydraulic fracture of the ductile rockusing a pore pressure cohesive zone model (CZM).
 9. The method of claim1, further comprising determining a flow of hydrocarbon in a fracturedductile rock.
 10. The method of claim 7, further comprising adjusting aparameter in a fracturing process based, at least in part, on thepredicted fracture in the ductile rock.
 11. The method of claim 10,wherein the parameter is a mass flow rate of a fluid, a pressure of afluid, or a combination thereof.
 12. A method for modeling hydraulicfractures in ductile rocks, comprising: generating a constitutive modelfor a quasi-brittle material that incorporates unified creep-plasticity;performing a finite element analysis using the constitutive model;determining an effective fracture property for the quasi-brittlematerial; and predicting hydraulic fractures in the quasi-brittlematerial.
 13. The method of claim 12, further comprising adjusting ahydraulic fracturing parameter based, at least in part, on the fracturespredicted in the ductile rock.
 14. The method of claim 12, wherein aratio between a length of a fracture process zone and a length of anon-linear fracture zone at a tip of a fracture is between about 1:1 and1:1.5.
 15. A system for modeling fractures in a reservoir, comprising: aprocessor; a storage medium comprising a data structure comprising amodel that incorporates unified creep-plasticity (UCP) into aconstitutive model for a ductile rock; and a computer readable mediumcomprising code configured to direct the processor to perform a finiteelement analysis of the ductile rock using the model.
 16. The system ofclaim 15, wherein the processor comprises a multi-processor clustercomputer system.
 17. The system of claim 15, wherein the computerreadable medium comprises code configured to direct the processor to runa finite element analysis for predicting a J-integral value and a stressfield incorporating plasticity effect using the model.
 18. The system ofclaim 15, wherein the computer readable medium comprises code configuredto direct the processor to determine the fracture properties of theductile rock.
 19. The system of claim 15, wherein the computer readablemedium comprises code configured to direct the processor to predict ahydraulic fracture in the ductile rock.
 20. The system of claim 15,wherein the computer readable medium comprises code configured to directthe processor to determine a flow of hydrocarbon in a fractured ductilerock.
 21. A non-transitory, computer-readable medium comprising a datastructure representing a cohesive fracture-mechanics model that combinesa constitutive model with a unified creep-plasticity model.
 22. Thenon-transitory, computer-readable medium of claim 21, further comprisingcode configured to direct a processor to: perform a finite elementanalysis using the constitutive model; determine an effective fractureproperty for a ductile rock; and predict fractures in the ductile rock.23. A method for harvesting hydrocarbons from a ductile rock formation,comprising: generating a model that incorporates unifiedcreep-plasticity (UCP) into a constitutive model for a ductile rock;performing a finite element analysis using the model; predictingfractures in a ductile rock formation based, at least in part, upon aresult from the finite element analysis; fracturing the ductile rockformation using parameters based, at least in part, on the prediction ofthe fractures; and harvesting hydrocarbons from the fractured rock.